The present invention relates to image processing systems and methods, and more particularly to image registration systems that combine two or more images into a composite image; in particular the fusion of anatomical manifold-based knowledge with volume imagery via large deformation mapping which supports both kinds of information simultaneously, as well as individually, and which can be implemented on a rapid convolution FFT based computer system.
Image registration involves combining two or more images, or selected points from the images, to produce a composite image containing data from each of the registered images. During registration, a transformation is computed that maps related points among the combined images so that points defining the same structure in each of the combined images are correlated in the composite image.
Currently, practitioners follow two different registration techniques. The first requires that an individual with expertise in the structure of the object represented in the images label a set of landmarks in each of the images that are to be registered. For example, when registering two MRI images of different axial slices of a human head, a physician may label points, or a contour surrounding these points, corresponding to the cerebellum in two images. The two images are then registered by relying on a known relationship among the landmarks in the two brain images. The mathematics underlying this registration process is known as small deformation multi-target registration.
In the previous example of two brain images being registered, using a purely operator-driven approach, a set of N landmarks identified by the physician, represented by xi, where i=1 . . . N, are defined within the two brain coordinate systems. A mapping relationship, mapping the N points selected in one image to the corresponding N points in the other image, is defined by the equation u(xi)=ki, where i=1 . . . N. Each of the coefficients, ki, is assumed known.
The mapping relationship u(x) is extended from the set of N landmark points to the continuum using a linear quadratic form regularization optimization of the equation:                     u        =                  arg          ⁢                      xe2x80x83                    ⁢                                    min              u                        ⁢                          xe2x80x83                        ⁢                          ∫                                                "LeftDoubleBracketingBar"                  Lu                  "RightDoubleBracketingBar"                                2                                                                        (        1        )            
subject to the boundary constraints u(xi)=kixe2x80x2. The operator L is a linear differential operator. This linear optimization problem has a closed form solution. Selecting L=xcex1∇2+xcex2∇(∇xc2x7) gives rise to small deformation elasticity.
For a description of small deformation elasticity see S. Timoshenko, Theory of Elasticity, McGraw-Hill, 1934 (hereinafter referred to as Timoshenko) and R. L. Bisplinghoff, J. W. Marr, and T. H. H. Pian, Statistics of Deformable Solids, Dover Publications, Inc., 1965 (hereinafter referred to as Bisplinghoff). Selecting L=∇2 gives rise to a membrane or Laplacian model. Others have used this operator in their work, see e.g., Amit, U. Grenander, and M. Piccioni, xe2x80x9cStructural image restoration through deformable templates,xe2x80x9d J. American Statistical Association. 86(414):376-387, June 1991, (hereinafter referred to as Amit) and R. Szeliski, Bayesian Modeling of Uncertainty in Low-Level Vision, Kluwer Academic Publisher, Boston, 1989 (hereinafter referred to as Szeliski) (also describing a bi-harmonic approach). Selecting L=∇4 gives a spline or biharmonic registration method. For examples of applications using this operator see Grace Wahba, xe2x80x9cSpline Models for Observational Data,xe2x80x9d Regional Conference Series in Applied Mathematics. SIAM, 1990, (hereinafter referred to as Whaba) and F. L. Bookstein, The Measurement of Biological Shape and Shape Change, volume 24, Springer-Verlag: Lecture Notes in Biomathematics, New York, 1978 (hereinafter referred to as Bookstein).
The second currently-practiced technique for image registration uses the mathematics of small deformation multi-target registration and is purely image data driven. Here, volume based imagery is generated of the two targets from which a coordinate system transformation is constructed. Using this approach, a distance measure, represented by the expression D(u), represents the distance between a template T(x) and a target image S(x). The optimization equation guiding the registration of the two images using a distance measure is:                     u        =                              arg            ⁢                          xe2x80x83                        ⁢                                          min                u                            ⁢                              xe2x80x83                            ⁢                              ∫                                                      "LeftDoubleBracketingBar"                    Lu                    "RightDoubleBracketingBar"                                    2                                                              +                      D            ⁡                          (              u              )                                                          (        2        )            
The distance measure D(u) measuring the disparity between imagery has various forms, e.g., the Gaussian squared error distance ∫|T(h(x))xe2x88x92S(x)|2dx, a correlation distance, or a Kullback Liebler distance. Registration of the two images requires finding a mapping that minimizes this distance.
Other fusion approaches involve small deformation mapping coordinates x xcex5 xcexa9 of one set of imagery to a second set of imagery. Other techniques include the mapping of predefined landmarks and imagery, both taken separately such as in the work of Bookstein, or fused as covered via the approach developed by Miller-Joshi-Christensen-Grenander the ""212 patent described in U.S. Pat. No. 6,009,212 (hereinafter xe2x80x9cthe ""212 patentxe2x80x9d) herein incorporated by reference mapping coordinates x xcex5 xcexa9 of one target to a second target. The existing state of the art for small deformation matching can be stated as follows:
Small Deformation Matching: Construct h(x)=xxe2x88x92u(x) according to the minimization of the distance D(u) between the template and target imagery subject to the smoothness penalty defined by the linear differential operator L:                               h          ⁡                      (            ·            )                          =                                            ·                              -                                  u                  ⁡                                      (                    ·                    )                                                                        ⁢                          xe2x80x83                        ⁢            where            ⁢                          xe2x80x83                        ⁢                          u              ^                                =                                    arg              ⁢                              xe2x80x83                            ⁢                                                min                  u                                ⁢                                  xe2x80x83                                ⁢                                  ∫                                                            "LeftDoubleBracketingBar"                      Lu                      "RightDoubleBracketingBar"                                        2                                                                        +                                          ∑                                  i                  =                  1                                N                            ⁢                                                D                  ⁡                                      (                    u                    )                                                  .                                                                        (        3        )            
The distance measure changes depending upon whether landmarks or imagery are being matched.
1. Landmarks alone. The first approach is purely operator driven in which a set of point landmarks xi, i=1, . . . N are defined within the two brain coordinate systems by, for example, an anatomical expert, or automated system from which the mapping is assumed known: u(xi)=ki, i=1, . . . , N. The field u(x) specifying the mapping h is extended from the set of points {xi} identified in the target to the points {yi} measured with Gaussian error co-variances xcexa3i:                               u          ^                =                              arg            ⁢                          xe2x80x83                        ⁢                                          min                u                            ⁢                              xe2x80x83                            ⁢                              ∫                                                      "LeftDoubleBracketingBar"                    Lu                    "RightDoubleBracketingBar"                                    2                                                              +                                    ∑                              i                =                1                            N                        ⁢                                                            (                                                            y                      i                                        -                                          x                      i                                        -                                          u                      ⁡                                              (                                                  x                          i                                                )                                                                              )                                t                            ⁢                                                                    Σ                                          -                      1                                                        ⁡                                      (                                                                  y                        i                                            -                                              x                        i                                            -                                              u                        ⁡                                                  (                                                      x                            i                                                    )                                                                                      )                                                  T                                                                        (        4        )            
Here (xc2x7)xcfx84 denotes transpose of the 3xc3x971 real vectors, and L is a linear differential operator giving rise to small deformation elasticity (see Timoshenko and Bisplinghoff), the membrane of Laplacian model (see Amit and Szeliski), bi-harmonic (see Szeliski), and many of the spline methods (see Wahba and Bookstein). This is a linear optimization problem with closed form solution.
2. The second approach is purely volume image data driven, in which the volume based imagery is generated of the two targets from which the coordinate system transformation is constructed. A distance measure between the two images being registered I0, I1 is defined as D(u)=∫|I0(xxe2x88x92u(x))xe2x88x92I1(x)|2dx. The corresponding optimization is:                                           h            ⁡                          (              ·              )                                =                                    ·                              -                                  u                  ⁡                                      (                    ·                    )                                                                        ⁢                          xe2x80x83                        ⁢            where                          ⁢                  xe2x80x83                ⁢                  
                ⁢                  u          =                                    arg              ⁢                              xe2x80x83                            ⁢                                                min                  u                                ⁢                                  xe2x80x83                                ⁢                                  ∫                                                            "LeftDoubleBracketingBar"                      Lu                      "RightDoubleBracketingBar"                                        2                                                                        +                          ∫                                                                    "LeftBracketingBar"                                                                                            I                          0                                                ⁡                                                  (                                                      x                            -                                                          u                              ⁡                                                              (                                x                                )                                                                                                              )                                                                    -                                                                        I                          1                                                ⁡                                                  (                          x                          )                                                                                      "RightBracketingBar"                                    2                                ⁢                                                      ⅆ                    x                                    .                                                                                        (        5        )            
The data function D(u) measures the disparity between imagery and has various forms. Other distances are used besides the Gaussian squared error distance, including correlation distance, Kullback Liebler distance, and others.
3. The algorithm for the transformation of imagery I0 into imagery I1 has landmark and volume imagery fused in the small deformation setting as in the ""212 patent. Both sources of information are combined into the small deformation registration:                     u        =                              arg            ⁢                          xe2x80x83                        ⁢                                          min                u                            ⁢                              xe2x80x83                            ⁢                              ∫                                                      "LeftDoubleBracketingBar"                    Lu                    "RightDoubleBracketingBar"                                    2                                                              +                      D            ⁡                          (              u              )                                +                                    ∑                              i                =                1                            N                        ⁢                                                            "LeftDoubleBracketingBar"                                                            y                      i                                        -                                          x                      i                                        -                                          u                      ⁡                                              (                                                  x                          i                                                )                                                                              "RightDoubleBracketingBar"                                2                                            σ                ⁢                                  xe2x80x83                                ⁢                                  2                  1                                                                                        (        6        )            
Although small deformation methods provide geometrically meaningful deformations under conditions where the imagery being matched are small, linear, or affine changes from one image to the other. Small deformation mapping does not allow the automatic calculation of tangents, curvature, surface areas, and geometric properties of the imagery. To illustrate the mapping problem FIG. 9 shows an oval template image with several landmarks highlighted. FIG. 10 shows a target image that is greatly deformed from the template image. The target image is a largely deformed oval that has been twisted. FIG. 11 shows the results of image matching when the four corners are fixed, using small deformation methods based on static quadratic form regularization. These figures illustrate the distortion which occurs with small deformation linear mapping when used with landmark points which define a motion corresponding to large deformation. As can be seen in FIG. 11, landmarks defined in the template image often map to more than one corresponding point in the target image.
Large deformation mapping produces maps for image registration in which the goal is to find the one-to-one, onto, invertible, differentiable maps h (henceforth termed diffeomorphisms) from the coordinates x xcex5 xcexa9 of one target to a second target under the mapping
h:xxe2x86x92h(x)=xxe2x88x92u(x),xxcex5xcexa9xe2x80x83xe2x80x83(7) 
To accommodate very fine variations in anatomy the diffeomorphic transformations constructed are of high dimensions having, for example a dimension greater than 12 of the Affine transform up-to the order of the number of voxels in the volume. A transformation is diffeomorphic if the transformation from the template to the target is one-to-one, onto, invertible, and both the transformation and it""s inverse are differentiable. A transformation is said to be one-to-one if no two distinct points in the template are mapped to the same point in the target. A transformation is said to be onto if every point in the target is mapped from a point in the template. The importance of generating diffeomorphisms is that tangents, curvature, surface areas, and geometric properties of the imagery can be calculated automatically. FIG. 12 illustrates the image mapping illustrated in FIG. 11 using diffeomorphic transformation.
The present invention overcomes the limitations of the conventional techniques of image registration by providing a methodology which combines, or fuses, some aspects of techniques where an individual with expertise in the structure of the object represented in the images labels a set of landmarks in each image that are to be registered and techniques that use mathematics of small deformation multi-target registration, which is purely image data driven.
Additional features and advantages of the invention will be set forth in the description which follows, and in part, will be apparent from the description, or may be learned by practicing the invention. The embodiments and other advantages of the invention will be realized and obtained by the method and apparatus particularly pointed out in the written description and the claims hereof as well as in the appended drawings.
To achieve these and other advantages and in accordance with the purpose of the invention, as embodied and broadly described, an embodiment consistent with the present invention registers a template and a target image by defining a manifold landmark point in the template image, identifying a point in said target image corresponding to the defined manifold landmark point, and selecting a coordinate frame suitable for spherical geometries. According to an embodiment of the present invention, a large deformation transform is computed using the selected coordinate frame, a manifold landmark transformation operator, and at least one manifold landmark transformation boundary value. The large deformation transform relates the manifold landmark point in the template image to the corresponding point in the target image. The template and target image are registered using the large deformation transform in the selected coordinate frame. Another embodiment consistent with the present invention registers a template image with a target image by selecting a coordinate frame suitable for spherical geometries and registering the target and template image using a large deformation transform in the selected coordinate frame.